install.packages("stargazer")

library(stargazer)


### Linear regression for cuts in public spending 
### (treatment effects without covariates)

I5_P11_1 <- lm(I5_P11_1 ~ wariant, 
               data = PBS_data)

I5_P11_2 <- lm(I5_P11_2 ~ wariant, 
               data = PBS_data)

I5_P11_3 <- lm(I5_P11_3 ~ wariant, 
               data = PBS_data)

I5_P11_4 <- lm(I5_P11_4 ~ wariant, 
               data = PBS_data)

I5_P11_5 <- lm(I5_P11_5 ~ wariant, 
   data = PBS_data)

I5_P11_6 <- lm(I5_P11_6 ~ wariant, 
               data = PBS_data)

I5_P11_7 <- lm(I5_P11_7 ~ wariant, 
               data = PBS_data)

I5_P11_8 <- lm(I5_P11_8 ~ wariant, 
               data = PBS_data)

I5_P11_se <- list(sqrt(diag(sandwich::vcovHC(I5_P11_1, type = "HC1"))),
                 sqrt(diag(sandwich::vcovHC(I5_P11_2, type = "HC1"))),
                 sqrt(diag(sandwich::vcovHC(I5_P11_3, type = "HC1"))),
                 sqrt(diag(sandwich::vcovHC(I5_P11_4, type = "HC1"))),
                 sqrt(diag(sandwich::vcovHC(I5_P11_5, type = "HC1"))),
                 sqrt(diag(sandwich::vcovHC(I5_P11_6, type = "HC1"))), 
                 sqrt(diag(sandwich::vcovHC(I5_P11_7, type = "HC1"))),
                 sqrt(diag(sandwich::vcovHC(I5_P11_8, type = "HC1"))))

lmtest::coeftest(I5_P11_1, df = Inf)

I5_P11 <- stargazer(I5_P11_1, I5_P11_2, I5_P11_3, I5_P11_4,
                 I5_P11_5, I5_P11_6, I5_P11_7, I5_P11_8,
                  header = FALSE, 
                  model.numbers = F,
                  omit.table.layout = "n",
                  digits = 3, 
                  column.labels = c("(1)", "(2)", "(3)", "(4)",
                                    "(5)", "(6)", "(7)", "(8)"),
                  dep.var.caption  = "Dependent Variable: Preference for cuts",
                  dep.var.labels.include = FALSE,
                  se = I5_P11_se)

### Regression for cuts in public spending
### With covariates to control for potential bias in randomization


I5_P11_1_c <- lm(I5_P11_1 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
                 + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_2_c  <- lm(I5_P11_2 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
                  + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_3_c  <- lm(I5_P11_3 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
                  + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_4_c  <- lm(I5_P11_4 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
            + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_5_c  <- lm(I5_P11_5 ~ wariant + M3_rekod + M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
            + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_6_c  <- lm(I5_P11_6 ~ wariant + M3_rekod + M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
            + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_7_c  <- lm(I5_P11_7 ~ wariant + M3_rekod + M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
            + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

I5_P11_8_c  <- lm(I5_P11_8 ~ wariant + M3_rekod + M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
            + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
           data = PBS_data)

lmtest::coeftest(I5_P11_4_c, df = Inf)

### Number of observation used in regression

nobs(I5_P11_8_c)

### Robust standard errors for regression with covariates

I5_P11_c_se <- list(sqrt(diag(sandwich::vcovHC(I5_P11_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I5_P11_2_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I5_P11_3_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I5_P11_4_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I5_P11_5_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I5_P11_6_c, type = "HC1"))), 
                sqrt(diag(sandwich::vcovHC(I5_P11_7_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I5_P11_8_c, type = "HC1"))))

### Latex code for regression 

I5_P11_c <- stargazer(I5_P11_1_c, I5_P11_2_c, I5_P11_3_c, I5_P11_4_c,
                      I5_P11_5_c, I5_P11_6_c, I5_P11_7_c, I5_P11_8_c,
          header = FALSE, 
          model.numbers = F,
          omit.table.layout = "n",
          digits = 3, 
          column.labels = c("(1)", "(2)", "(3)", "(4)",
                            "(5)", "(6)", "(7)", "(8)"),
          dep.var.caption  = "Dependent Variable: Preference for cuts",
          dep.var.labels.include = FALSE,
          se = I5_P11_c_se)

### Latex code for both tables (covariates and without covariates)
        
panelr <- starpolishr::star_panel(I5r, I5rc,
           panel.names = c("Without controls", "With controls"))
         
print(panelr)

### Regression for hikes in public spending
### Without covariates

I6_P12_1 <- lm(I6_P12_1 ~ wariant, 
           data = PBS_data)

I6_P12_2 <- lm(I6_P12_2 ~ wariant, 
           data = PBS_data)

I6_P12_3 <- lm(I6_P12_3 ~ wariant, 
           data = PBS_data)

I6_P12_4 <- lm(I6_P12_4 ~ wariant, 
           data = PBS_data)

I6_P12_5 <- lm(I6_P12_5 ~ wariant, 
           data = PBS_data)

I6_P12_6 <- lm(I6_P12_6 ~ wariant, 
           data = PBS_data)

I6_P12_7 <- lm(I6_P12_7 ~ wariant, 
           data = PBS_data)

I6_P12_8 <- lm(I6_P12_8 ~ wariant, 
           data = PBS_data)

I6_P12_8_se <- list(sqrt(diag(sandwich::vcovHC(I6_P12_1, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_2, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_3, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_4, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_5, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_6, type = "HC1"))), 
                sqrt(diag(sandwich::vcovHC(I6_P12_7, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_8, type = "HC1"))))

I6_P12 <- stargazer(I6_P12_1, I6_P12_2, I6_P12_3, I6_P12_4,
                 I6_P12_5, I6_P12_6, I6_P12_7, I6_P12_8,
                 header = FALSE, 
                 model.numbers = F,
                 omit.table.layout = "n",
                 digits = 3, 
                 column.labels = c("(1)", "(2)", "(3)", "(4)",
                                   "(5)", "(6)", "(7)", "(8)"),
                 dep.var.caption  = "Dependent Variable: Preference for cuts",
                 dep.var.labels.include = FALSE,
                 se = I6_P12_8_se)

### Regression for incerase in public spending
### With covariates to control for potential bias in randomization


I6_P12_1_c <- lm(I6_P12_1 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M11_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_2_c <- lm(I6_P12_2 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_3_c <- lm(I6_P12_3 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_4_c <- lm(I6_P12_4 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_5_c <- lm(I6_P12_5 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_6_c <- lm(I6_P12_6 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_7_c <- lm(I6_P12_7 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)

I6_P12_8_c <- lm(I6_P12_8 ~ wariant + M3_rekod +  M1 + M2_rekod + M3_rekod + M5 + M6 + M13 + M15 + M16_p + M17 
               + M_20 + M11_1_b + M12_2_b + M11_3_b + M11_5_b, 
            data = PBS_data)


lmtest::coeftest(I6_P12_1_c, df = Inf)

### Number of observation used in regression

nobs(I61rc)

### Robust standard errors for regression with covariates

I6_P12_se_c <- list(sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))), 
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))),
                sqrt(diag(sandwich::vcovHC(I6_P12_1_c, type = "HC1"))))

### Latex code for regression 

I6rc <- stargazer(I6_P12_1_c, I6_P12_2_c, I6_P12_3_c, I6_P12_4_c,
                  I6_P12_5_c, I6_P12_6_c, I6_P12_7_c, I6_P12_8_c,
                  header = FALSE, 
                  model.numbers = F,
                  omit.table.layout = "n",
                  digits = 3, 
                  column.labels = c("(1)", "(2)", "(3)", "(4)",
                                    "(5)", "(6)", "(7)", "(8)"),
                  dep.var.caption  = "Dependent Variable: Preference for cuts",
                  dep.var.labels.include = FALSE,
                  se = I6_P12_se_c)

